Random intercept model met 2 niveaus
Op deze pagina vind je een demonstratie van een statistische techniek aan de hand van een voorbeeld.
Meer informatie over hoe je deze pagina kan gebruiken vind je in deze handleiding.
De analyse gebeurt met behulp van R en RStudio. Een inleiding tot deze software vind je hier.
Het voorbeeld op deze pagina is afkomstig van een studie van Bogaert, Merchie, Van Ammel & Van Keer (2025). Er zijn enkele wijzigingen aangebracht om didactische redenen. De studie is ook uitgebreider dan wat op deze pagina wordt gedemonstreerd.
1 Het voorbeeld op deze pagina
Bogaert et al. (2025) onderzoeken de impact van een interventie op verschillende maten van leesvaardigheid van leerlingen. Hun studie is gemotiveerd door de vaststelling dat leesvaardigheden vaak onvoldoende ontwikkeld zijn bij leerlingen in de latere jaren van het basisonderwijs. Bovendien zijn er onbeantwoorde vragen over welke lesmethoden effectief zouden kunnen zijn bij het verbeteren van de leesvaardigheden van leerlingen. Studies die uitzoeken of interventies op klasniveau zouden kunnen werken zijn bijzonder schaars. De auteurs hadden de intentie om die lacune in de bestaande wetenschappelijke literatuur te verhelpen.
De interventie in de studie van Bogaert et al. (2025) houdt in dat leerlingen uit het basisonderwijs gedurende tien weken een lessenpakket volgen dat specifiek gericht is op het verbeteren van verschillende uitkomsten:
- de mate waarin leerlingen een tekst begrijpen (“reading comprehension”)
- het gebruik van verschillende strategieën om een tekst te begrijpen (“reading strategy use”)
- verschillende vormen van leesmotivatie (“reading motivation”)
In de demonstratie op deze pagina focussen we enkel op de eerste van die uitkomstvariabelen: het begrijpend lezen.
Naast de leerlingen bij wie de interventie plaatsvond, bestudeerden
Bogaert et al. (2025) ook leerlingen bij wie geen aanpassingen aan het
lessenpakket zijn gebeurd. Zij hebben lessen gevolgd zoals gewoonlijk.
In de variabele Conditie wordt geregistreerd of een
leerling tot de experimentele groep (met interventie) of tot de
controlegroep (zonder interventie) behoort.1
Alle leerlingen legden twee testen af die peilen naar capaciteiten
met betrekking tot begrijpend lezen. De eerste test vond plaats voor de
interventie (variabele BL_pretest), de tweede erna
(variabele BL_posttest).2
Een belangrijke vraag die de auteurs zich stelden was of de
interventie (Conditie) gemiddeld genomen een effect heeft
op de testscore (BL_posttest), rekening houdend met
eventuele andere predictoren.
Op het eerste zicht lijkt dit een vraag die kan worden beantwoord aan
de hand van een klassiek
lineair
regressiemodel met een categorische predictor
(Conditie). In de volgende sectie zetten we kort uiteen
waarom dat niet het geval is en hoe het beter kan.
2 Waarom een multilevel model met random intercept?
Om te begrijpen waarom klassieke lineaire regressie hier niet geschikt is moeten we focussen op de manier waarop de data zijn verzameld. In Bogaert et al. (2025) werden willekeurig scholen geselecteerd en vervolgens klassen in die scholen en daarna werden leerlingen binnen die klassen getest. De data zijn met andere woorden hiërarchisch gestructureerd, in dit geval met drie niveaus: scholen, klassen en individuele leerlingen. Voor deze inleidende demonstratie zullen we de context wat vereenvoudigen naar een structuur met twee niveaus: we gaan ervan uit dat willekeurige klassen werden geselecteerd en dat vervolgens leerlingen binnen die klassen werden getest.
Als onderzoeker zou je kunnen vermoeden dat de uitkomst (testscores
BL_posttest) geclusterd is per klas. Clustering
betekent dat de scores binnen klassen meer op elkaar lijken dan wanneer
je de scores over alle klassen heen zou beschouwen. Op basis van de klas
waartoe een leerling behoort zou je dus al iets kunnen zeggen over de
score die die leerling behaalt op de test (BL_posttest). Er
zou met andere woorden een “klaseffect” zijn op de testscore.3
De mogelijke aanwezigheid van een klaseffect kan relevant zijn bij
het beantwoorden van de onderzoeksvraag. De interesse van de
onderzoekers gaat uit naar het effect van Conditie op
BL_posttest. Als we een zuiver beeld willen van dit effect,
dan moeten we controleren voor het klaseffect. De doelstelling is
dezelfde als die van statistische controle bij lineaire regressie. Er
zijn echter goede redenen om het niet op die klassieke manier aan te
pakken.
De belangrijkste reden is dat de niveaus van de categorische
variabele Klas op een willekeurige wijze aanwezig zijn in
de data. Dat is een direct gevolg van de manier waarop de steekproef tot
stand is gekomen: de klassen zijn willekeurig getrokken uit alle
mogelijke klassen. De waarden van Klas die in een
steekproef voorkomen zijn een willekeurige greep uit alle mogelijke
waarden. In jouw steekproef zal misschien “het 6e studiejaar van
Stedelijke Basisschool De Brug in Roeselare” wel aanwezig zijn en “het
5e studiejaar van Vrije Basisschool De Klimop in Sint-Gillis-Waas” niet.
In een andere steekproef zou het omgekeerd kunnen zijn.
Dat is natuurlijk heel verschillend van wat je gewoon bent bij
categorische variabelen. Bijvoorbeeld, bij een variabele als
Opleidingsniveau met 4 niveaus zal je niet een willekeurige
trekking van die niveaus nemen, om vervolgens binnen elk getrokken
niveau willekeurig individuen te trekken.
Doordat de niveaus van Klas willekeurig getrokken zijn,
zijn ook de effecten die geassocieerd zijn met die niveaus willekeurig.
Het zijn met andere woorden random effecten. Als je in een
dergelijke context een klassiek lineair regressiemodel zou fitten, dan
kan je je conclusies maar in beperkte mate generaliseren. Je conclusies
gelden voor alle leerlingen binnen de specifieke getrokken
klassen.4 Een multilevel model daarentegen laat toe
om je conclusies wél te generaliseren naar alle leerlingen in alle
klassen.
Een bijkomende reden om voor een multilevel model te kiezen is
efficiënt gebruik van de data. De dataset kan heel veel klassen
bevatten. Daardoor zou een lineair regressiemodel heel veel parameters
bevatten waar we niet echt in geïnteresseerd zijn, maar die wel moeten
worden geschat. Bijvoorbeeld, het effect van in klas “het 6e studiejaar
van Stedelijke Basisschool De Brug in Roeselare” te zitten in
vergelijking met het referentieniveau. Analoog zou je het effect moeten
schatten voor alle klassen in je dataset. Dat zouden we moeten doen om
de statistische controle te kunnen uitvoeren, niet omdat al die effecten
ons inhoudelijk interesseren. Dat is geen efficiënt gebruik van de data,
die beter zouden worden gespendeerd aan het verbeteren van de power van
de toets waar we wel in geïnteresseerd zijn, in dit geval de toets van
het effect van Conditie.
In deze demonstratie komen enkel random intercepten aan bod. Daarbij
verschilt dus het niveau van de uitkomst BL_posttest
(random) per klas, ongeacht de Conditie. Merk op dat er ook
“random slope” modellen bestaan. Daarbij kan niet alleen het intercept,
maar ook de sterkte van het verband (en dus de helling van de rechte)
variëren per klas. Dergelijke modellen worden op deze pagina niet
gedemonstreerd. Een goede maar erg beknopte Engelstalige demonstratie
vind je
hier.
Terminologie
Multilevel modellen staan ook bekend als (linear) mixed models, mixed effects models en hiërarchische lineaire modellen.
3 Data & packages
Voor een multilevel analyse in R gebruikt men meestal twee packages:
lme4 en lmerTest. Indien deze nog niet
geïnstalleerd zijn op je computer dan moet je dit eerst eenmalig
doen.
install.packages(c('lme4', 'lmerTest'))
Als de packages geïnstalleerd zijn laad je ze voor gebruik. Het is in
dit geval belangrijk dat je de volgorde van het laden respecteert: eerst
lme4, dan lmerTest.5
library(lme4)
library(lmerTest)
Verderop zullen we voor enkele heel specifieke doeleinden nog een aantal packages gebruiken. We verduidelijken dit wanneer het relevant is.
De data voor deze demonstratie komen van Bogaert et al. (2025). Je kan ze inladen met de volgende code. De namen van de variabelen zijn vertaald en licht aangepast voor de duidelijkheid.
lezen <- read.csv('https://statlas.ugent.be/datasets/mlm.csv')
Het is een goed idee om categorische variabelen meteen om te zetten
naar het datatype factor. Indien gewenst kan je ook
labels toevoegen. Uiteraard is het belangrijk om hier geen
fouten te maken door de volgorde om te keren.
lezen$School <- factor(lezen$School)
lezen$Klas <- factor(lezen$Klas)
lezen$Conditie <- factor(lezen$Conditie, labels=c('Controle', 'Experimenteel'))
lezen$Niveau <- factor(lezen$Niveau, labels=c('Goed', 'Middelmatig', 'Zwak'))
lezen$Dyslexie <- factor(lezen$Dyslexie, labels=c('Nee', 'Ja'))
4 Is de uitkomst geclusterd?
Het kan interessant zijn om te weten in welke mate de uitkomstvariabele geclusterd is op een hoger niveau.
In Bogaert et al. (2025) ligt de focus op de uitkomstvariabele
BL_posttest. Deze variabele is gemeten op het niveau van de
individuele leerlingen (niveau 1).
Leerlingen maken deel uit van klassen (niveau 2). Als onderzoeker kan je vermoeden dat klassen van elkaar verschillen in termen van scores op de posttest. Anders gezegd, de posttest-scores binnen klassen lijken mogelijk wat meer op elkaar dan wanneer je alle leerlingen over alle klassen heen zou beschouwen.
We vragen ons af of dat vermoeden klopt. Om dat te beantwoorden
bouwen we een nulmodel of een intercept-only model. Dat dient enkel om
de variantie in BL_posttest op te splitsen in variantie op
niveau 1 en variantie op niveau 2. De functie lmer() laat
toe om een nulmodel te bouwen. Bij het argument formula
specifiëren we een model zonder echte predictoren. Door
(1|Klas) te schrijven laten we wel een afzonderlijk
intercept toe per klas. Met summary() vragen we de
belangrijkste informatie over dit model op.
model0 <- lmer(formula = BL_posttest ~ (1|Klas), data=lezen)
summary(model0) Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: BL_posttest ~ (1 | Klas)
Data: lezen
REML criterion at convergence: 1579.3
Scaled residuals:
Min 1Q Median 3Q Max
-2.5693 -0.6433 -0.1237 0.4940 7.5250
Random effects:
Groups Name Variance Std.Dev.
Klas (Intercept) 0.1928 0.4391
Residual 1.6473 1.2835
Number of obs: 464, groups: Klas, 27
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.02125 0.10464 26.50000 0.203 0.841In de output van summary(model0), onder
Random effects vinden we de opgesplitste variantie terug.
De variantie op niveau 1 (naast Residual) is duidelijk het
grootst, maar er is toch ook variantie op het niveau van de klassen.
Een visualisatie van de data opgesplitst per klas kan helpen om een idee te krijgen van de verhouding tussen deze varianties.
boxplot(lezen$BL_posttest ~ lezen$Klas)
Elke boxplot stelt de data in een van de 27 klassen voor. We komen tot een gelijkaardige conclusie: er lijkt toch variantie te bestaan tussen de klassen, met andere woorden variantie op niveau 2.
Om de verhouding tussen de varianties op verschillende niveaus te kwantificeren kunnen we de intraclass correlatie (ICC) berekenen.
4.1 Intraclass correlatie (ICC)
De intraclass correlatie (ICC) drukt uit hoe sterk de clustering op
niveau 2 is. De ICC geeft de proportie weer van de variantie op niveau 2
(\(\sigma^2_{klas}\)) ten opzichte van
de totale variantie (\(\sigma^2_{klas} +
\sigma^2_{res}\)). De ICC kan waarden aannemen van 0 tot 1.6 De nodige
gegevens om de ICC te berekenen zijn te vinden in de output van
summary(model0).
\[ICC =
\frac{\sigma^2_{klas}}{\sigma^2_{klas} + \sigma^2_{res}}\]
In onze steekproef wordt dat
\[\frac{0.1928}{0.1928 + 1.6473} = 0.1049823 \approx 0.105\]
In dit geval is de ICC ongeveer gelijk aan \(0.105\) of \(10.5\%\).
Het package performance biedt een functie
icc aan die de berekening voor ons kan maken. Installeer
het package met install.packages() indien je dat nog niet
hebt gedaan en laad het met library().
install.packages('performance')
library(performance)
icc(model0)# Intraclass Correlation Coefficient
Adjusted ICC: 0.105
Unadjusted ICC: 0.105Als we een nulmodel aan de functie icc() geven zal er
geen verschil zijn tussen de Adjusted en de
Unadjusted varianten. Met het commando ?icc
kan je de helppagina opvragen waar je een uitleg vindt over het verschil
tussen beide.
4.2 Random intercepten visualiseren
Een makkelijke manier om de random intercepten te visualiseren is met
de functie dotplot() van het package lattice.
Met de functie ranef() kan je de random intercepten van een
model (hier model0) opvragen. Die geef je aan de functie
dotplot(). Je ziet in de output voor elke klas het random
intercept.
library(lattice)
dotplot(ranef(model0))
Multilevel modellen maken de assumptie dat random intercepten een
normale verdeling volgen. Met de functies qqnorm() en
qqline() is het mogelijk om na te gaan of aan deze
assumptie is voldaan. Hieronder zie je dit toegepast voor het nulmodel
model0 en voor clustervariabele Klas (de enige
clustervariabele in dit voorbeeld).
qqnorm(ranef(model0)$Klas[,1])
qqline(ranef(model0)$Klas[,1])
De punten wijken niet heel veel af van de rechte. Bijgevolg hebben we geen overtuigende reden om te denken dat de assumptie van normaliteit geschonden is.
5 Effecten schatten en toetsen
Bogaert et al. (2025) zijn geïnteresseerd in het effect van een
interventie (categorische variabele Conditie) op de
testscore van leerlingen na de interventie BL_posttest. Een
bijkomende controlevariabele in hun model is de score van leerlingen op
de test vóór de interventie (continue variabele
BL_pretest).
Algemeen gesteld kan je gelijk welke combinatie van categorische en continue predictoren - ook met interactie-effecten - in een multilevel model opnemen, net als bij klassieke lineaire regressie.
We analyseren hieronder eerst het effect van BL_pretest
op BL_posttest en daarna voegen we de predictor
Conditie toe aan het model. Dit laat ons toe om zowel de
analyse van een continue predictor als van een categorische predictor te
demonstreren. Daarna tonen we hoe je een interactie-effect kan
toevoegen.
5.1 Continue predictor
BL_pretest
Eerst bouwen we het model met enkel de pretestscores
BL_pretest als predictor, naast het random intercept per
Klas. Het effect van deze pretestscores is wetenschappelijk
niet bijzonder interessant. Het is ook niet waarin Bogaert et al. (2025)
direct geïnteresseerd zijn. We staan er enkel bij stil om te
demonstreren hoe je het effect van een continue predictor kan
toetsen.
Specificatie in R
De specificatie van het model met predictor BL_pretest
is gelijkaardig aan de specificatie van het nulmodel. Het enige verschil
is de toevoeging van BL_pretest bij het argument
formula.
model.prtst <- lmer(formula = BL_posttest ~ BL_pretest + (1|Klas), data=lezen)
In een diagram ziet dit model er zo uit (klik op de afbeelding om te vergroten):
×
Parameters schatten
De belangrijkste informatie over de parameterschattingen bekomen we
met de functie summary().
summary(model.prtst)Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: BL_posttest ~ BL_pretest + (1 | Klas)
Data: lezen
REML criterion at convergence: 1442.6
Scaled residuals:
Min 1Q Median 3Q Max
-2.3733 -0.6071 -0.0867 0.4709 7.7703
Random effects:
Groups Name Variance Std.Dev.
Klas (Intercept) 0.02737 0.1654
Residual 1.34285 1.1588
Number of obs: 456, groups: Klas, 27
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.02334 0.06347 22.12260 0.368 0.717
BL_pretest 0.51553 0.04245 402.82000 12.144 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
BL_pretest -0.006
In deze output is het de coëfficiënt onder Fixed effects
bij BL_pretest die ons het meest interesseert.
We zien daar 0.51553 staan. Dat is de puntschatting voor het effect
van BL_pretest op BL_posttest, controlerend
voor het klaseffect. Hier schatten we dat de score op
BL_posttest, gemiddeld over alle leerlingen van alle
klassen heen, 0.51553 punten hoger ligt wanneer de score op
BL_pretest met 1 eenheid toeneemt, controlerend voor het
klaseffect.
Het is een goed idee om naast puntschattingen ook
betrouwbaarheidsintervallen te rapporteren. Dat kan met de functie
confint(). Standaard krijg je een
95%-betrouwbaarheidsinterval.
confint(model.prtst)Computing profile confidence intervals ... 2.5 % 97.5 %
.sig01 0.0000000 0.3306104
.sigma 1.0843316 1.2399720
(Intercept) -0.1033102 0.1497472
BL_pretest 0.4301067 0.6052039
Terzijde, de output van summary(model.prtst) bevat ook
een intercept bij de Fixed effects. Het gaat hier over het
intercept over alle klassen heen. Voor een leerling met
BL_pretest gelijk aan 0 verwachten we een score van
0.02334 op de BL_posttest.
Toetsen
De relevante nulhypothese en alternatieve hypothese zijn de volgende.
-
De nulhypothese stelt dat de populatieparameter \(\beta_{pretest}\) gelijk is aan 0. In dat
geval is er geen effect van
BL_pretest. -
De alternatieve hypothese stelt dat de populatieparameter \(\beta_{pretest}\) verschillend is van 0. In
dat geval is er wel een effect van
BL_pretest.
Op dezelfde rij in de output van summary() (zie
hierboven) vinden we ook de nodige informatie om het effect van
BL_pretest te toetsen. Deze werkwijze zal echter niet in
elke situatie werken. Met name bij een categorische predictor en bij
modellen met een interactie-effect zal je met summary()
vaak niet de gewenste toets kunnen uitvoeren. Daarom demonstreren we
hier meteen een werkwijze die wel algemeen toepasbaar is.
Die werkwijze bestaan erin dat we:
- overgaan tot effectcodering van alle categorische predictoren
-
een modelvergelijking uitvoeren met de functie
Anova()uit het packagecar -
het argument
type=3meegeven aanAnova()
Meer uitleg over effectcodering en type III toetsen kan je vinden in deze syllabus vanaf sectie 5.3. De uitleg gebeurt in de context van klassieke lineaire regressie (dus niet multilevel modellen) maar de principes zijn identiek.
In model.prtst zijn geen categorische predictoren, dus
die stap kunnen we overslaan. Het argument type=3 speelt
enkel een rol wanneer er een interactie-effect in het model zit, maar
voor de volledigheid voegen we het hier toch toe.
library(car)
Anova(model.prtst, type=3) # Let op de hoofdletter "A"Analysis of Deviance Table (Type III Wald chisquare tests)
Response: BL_posttest
Chisq Df Pr(>Chisq)
(Intercept) 0.1352 1 0.7131
BL_pretest 147.4707 1 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Aangezien de p-waarde ongeveer gelijk is aan 0 kunnen we de nulhypothese verwerpen op het 5% significantieniveau.
Voor een vollediger voorbeeld van deze werkwijze, zie verderop.
Visualisatie
Het model dat we zonet hebben gefit met lmer() kan
visueel worden weergegeven als parallelle rechten door een puntenwolk.7 Er is 1
rechte per klas. De plot is wat ingezoomd om de verschillende rechten
iets duidelijker weer te geven.
- Elke klas heeft een eigen rechte omdat elke klas een eigen random intercept heeft.
-
De rechten zijn parallel omdat het effect van
BL_pretestin elke klas gelijk is. Het is met andere woorden louter een fixed effect in dit model.8
Bovenstaande visualisatie is omslachtig om te maken en leggen we daarom hier niet verder uit.
Een visualisatie enkel van het fixed effect van
BL_pretest is makkelijker, dankzij het package
effects. Je geeft de predictor en het model aan de functie
effect() en geeft dit op zijn beurt aan de functie
plot().
install.packages('effects') # eenmalig
library(effects)
plot(effect('BL_pretest', model.prtst)) 
5.2 Categorische
predictor Conditie
De centrale focus in Bogaert et al. (2025) was een hypothese over het
effect van de predictor Conditie op
BL_posttest.
Specificatie in R
We willen een model met de predictor Conditie, samen met
BL_pretest en het random intercept. Dat kan heel eenvoudig
door een term Conditie toe te voegen aan het argument
formula.
model.cndt <- lmer(formula = BL_posttest ~ BL_pretest + Conditie + (1|Klas), data=lezen)
In een diagram ziet dit model er zo uit (klik op de afbeelding om te vergroten):
×
Merk op dat de predictor Conditie in het diagram
gesitueerd is op niveau 2. De interventie zelf bestaat uit een
lessenpakket dat klassikaal wordt gebracht aan de leerlingen. Daardoor
behoren alle leerlingen van dezelfde klas onvermijdelijk tot dezelfde
conditie. De predictor Conditie is dus eerder een variabele
op klasniveau dan op het individuele niveau. Dat hoeft geen probleem te
zijn. Het is perfect mogelijk om predictoren van verschillende niveaus
te combineren in een multilevel model.
Parameters schatten
We vragen de belangrijkste informatie met betrekking tot de
parameterschattingen op met de functie summary(). De
predictor Conditie is een binaire categorische variabele.
Daardoor is er slechts 1 parameter die voor het effect van
Conditie codeert. De schatting van die ene parameter is te
vinden onder Fixed effects bij
ConditieExperimenteel.
summary(model.cndt)Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: BL_posttest ~ BL_pretest + Conditie + (1 | Klas)
Data: lezen
REML criterion at convergence: 1444.8
Scaled residuals:
Min 1Q Median 3Q Max
-2.3758 -0.6082 -0.0870 0.4655 7.7631
Random effects:
Groups Name Variance Std.Dev.
Klas (Intercept) 0.03204 0.179
Residual 1.34264 1.159
Number of obs: 456, groups: Klas, 27
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.02935 0.09190 20.81138 0.319 0.753
BL_pretest 0.51366 0.04258 404.54484 12.065 <2e-16 ***
ConditieExperimenteel -0.01206 0.12980 21.32136 -0.093 0.927
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) BL_prt
BL_pretest 0.009
CndtExprmnt -0.708 -0.018
Hier wordt geschat dat het gemiddelde effect, over alle leerlingen van alle klassen heen, van in de experimentele groep te zitten in vergelijking met de controlegroep (= het referentieniveau), controlerend voor de andere predictoren, gelijk is aan -0.01206.
Een leerling in de experimentele conditie heeft een score op
BL_posttest die, gemiddeld over alle leerlingen van alle
klassen heen, 0.01206 punten lager ligt dan een leerling in de
controleconditie, controlerend voor het klaseffect en voor
BL_pretest.
Het is een goed idee om naast puntschattingen ook
betrouwbaarheidsintervallen voor de parameters te rapporteren. Dat kan
met de functie confint(). Standaard krijg je een
95%-betrouwbaarheidsinterval.
confint(model.cndt)Computing profile confidence intervals ... 2.5 % 97.5 %
.sig01 0.0000000 0.3302957
.sigma 1.0843798 1.2400675
(Intercept) -0.1506928 0.2077439
BL_pretest 0.4302338 0.6055554
ConditieExperimenteel -0.2637086 0.2424676
Toetsen
De relevante nulhypothese en alternatieve hypothese bij het toetsen van het effect van een categorische predictor zijn de volgende.
-
De nulhypothese stelt dat alle populatieparameters (“bèta’s”) die
coderen voor de categorische variabele
Conditiegelijk zijn aan 0. Dat betekent dat er geen effect is vanConditie. -
De alternatieve hypothese stelt dat minstens een van deze
populatieparameters verschillend is van 0. Als dat klopt is er wel een
effect van
Conditie.
In dit specifieke geval zouden we de toets kunnen uitvoeren op basis
van de output van summary().9 Omdat dit in veel
situaties niet mogelijk is, zullen we hier focussen op een algemeen
toepasbare werkwijze:
- overgaan tot effectcodering van alle categorische predictoren
-
een modelvergelijking uitvoeren met de functie
Anova()uit het packagecar -
het argument
type=3meegeven aanAnova()
Meer uitleg over effectcodering en type III toetsen kan je vinden in deze syllabus vanaf sectie 5.3. De uitleg gebeurt in de context van klassieke lineaire regressie (dus niet multilevel modellen) maar de principes zijn identiek.
Voor de eerste stap moeten we het model opnieuw fitten, waarbij we
een extra argument contrasts meegeven.
model.cndt.effect <- lmer(formula = BL_posttest ~ BL_pretest + Conditie + (1|Klas),
contrasts=list(Conditie=contr.sum),
data=lezen)
library(car)
Anova(model.cndt.effect, type=3)Analysis of Deviance Table (Type III Wald chisquare tests)
Response: BL_posttest
Chisq Df Pr(>Chisq)
(Intercept) 0.1291 1 0.7194
BL_pretest 145.5551 1 <2e-16 ***
Conditie 0.0086 1 0.9260
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Het argument type=3 speelt enkel een rol wanneer er een
interactie-effect in het model aanwezig is, maar voor de volledigheid
voegen we het hier toch toe.
De p-waarde is ongeveer gelijk aan 0.926. Dit is groter dan 0.05 waardoor we de nulhypothese niet kunnen verwerpen op het 5% significantieniveau.
Visualisatie
Het (fixed) effect van Conditie kan eenvoudig
gevisualiseerd worden met het package effects. Je geeft de
predictor en het model aan de functie effect() en geeft dit
op zijn beurt aan de functie plot().
install.packages('effects') # eenmalig
library(effects)
plot(effect('Conditie', model.cndt))
5.3 Een interactie-effect
Bogaert et al. (2025) vermoeden dat het effect van de interventie
(Conditie) op de posttestscore (BL_posttest)
verschillend zou kunnen zijn voor leerlingen met dyslexie in
vergelijking met leerlingen zonder dyslexie. Daarom beslissen zij om een
interactie tussen de variabelen Conditie en
Dyslexie aan het model toe te voegen.
Specificatie in R
Met een asterisk * kan je een interactie-effect
toevoegen aan het model. De asterisk zorgt ervoor dat, naast het
interactie-effect zelf, de variabelen Conditie en
Dyslexie allebei ook als hoofdeffect in het model worden
opgenomen.10 Dat is wenselijk.
model.ia <- lmer(formula = BL_posttest ~ BL_pretest + Conditie*Dyslexie + (1|Klas), data=lezen)
In een diagram ziet dit model er zo uit (klik op de afbeelding om te vergroten):
×
Parameters schatten
De functie summary() geeft ons de belangrijkste
informatie over de parameterschattingen.
summary(model.ia)Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: BL_posttest ~ BL_pretest + Conditie * Dyslexie + (1 | Klas)
Data: lezen
REML criterion at convergence: 1364.1
Scaled residuals:
Min 1Q Median 3Q Max
-2.3667 -0.6204 -0.1125 0.4693 7.6658
Random effects:
Groups Name Variance Std.Dev.
Klas (Intercept) 0.03165 0.1779
Residual 1.37429 1.1723
Number of obs: 428, groups: Klas, 27
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.04952 0.09569 22.62317 0.517 0.610
BL_pretest 0.51244 0.04457 373.36755 11.497 <2e-16 ***
ConditieExperimenteel -0.01355 0.13571 23.79687 -0.100 0.921
DyslexieJa -0.29394 0.38402 422.60242 -0.765 0.444
ConditieExperimenteel:DyslexieJa 0.27432 0.57598 422.60653 0.476 0.634
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) BL_prt CndtEx DyslxJ
BL_pretest -0.008
CndtExprmnt -0.705 -0.016
DyslexieJa -0.189 0.017 0.133
CndtExpr:DJ 0.126 -0.061 -0.170 -0.668
Omdat zowel Conditie als Dyslexie slechts 2
waarden kunnen aannemen, is er 1 parameter die codeert voor het
interactie-effect. De schatting van deze parameter vinden we onder
Fixed effects, naast
ConditieExperimenteel:DyslexieJa.
De interpretatie van de parameter is dezelfde als bij een interactie-effect in een klassiek lineair regressiemodel. Het is het verschil in het effect van de ene variabele in de interactie op de uitkomst, wanneer de andere variabele in de interactie met 1 eenheid toeneemt en wanneer alle andere predictoren constant blijven.
In dit geval schatten we dat het effect van de interventie op
BL_posttest gemiddeld over alle leerlingen van alle klassen
heen, toeneemt met 0.27432 wanneer de dummyvariabele die codeert voor
Dyslexie met 1 eenheid toeneemt, terwijl
BL_pretest en Klas gelijk blijven. Wanneer de
dummy voor Dyslexie toeneemt met 1 eenheid betekent dit dat
Dyslexie de waarde Ja heeft in plaats van
Nee.11
Technisch gezien kan de schatting even goed geïnterpreteerd worden in
de andere richting, dus als het verschil in het effect van
Dyslexie op BL_posttest wanneer
Conditie met 1 eenheid toeneemt. Gegeven de hypothese van
de onderzoekers in dit voorbeeld is dat echter niet de relevante
interpretatie.
Het is een goed idee om naast puntschattingen ook
betrouwbaarheidsintervallen voor de parameters te rapporteren. Dat kan
met de functie confint().
confint(model.ia)Computing profile confidence intervals ... 2.5 % 97.5 %
.sig01 0.0000000 0.3331489
.sigma 1.0920458 1.2542892
(Intercept) -0.1379511 0.2344323
BL_pretest 0.4251270 0.6083951
ConditieExperimenteel -0.2755231 0.2525781
DyslexieJa -1.0442005 0.4552143
ConditieExperimenteel:DyslexieJa -0.8531328 1.3959948
Toetsen
De relevante nulhypothese en alternatieve hypothese bij het toetsen van een interactie-effect zijn de volgende.
-
De nulhypothese stelt dat de populatieparameter \(\beta_{Conditie \times Dyslexie}\) gelijk
is aan 0. In dat geval is er geen interactie-effect tussen
ConditieenDyslexie. -
De alternatieve hypothese stelt dat de populatieparameter \(\beta_{Conditie \times Dyslexie}\)
verschillend is van 0. In dat geval is er wel een interactie-effect
tussen
ConditieenDyslexie.
In dit specifieke geval zouden we de toets van het interactie-effect
kunnen uitvoeren op basis van de output van summary().12 Omdat
dit in veel situaties niet mogelijk is, zullen we hier focussen op een
algemeen toepasbare werkwijze:
- overgaan tot effectcodering van alle categorische predictoren
-
een modelvergelijking uitvoeren met de functie
Anova()uit het packagecar -
het argument
type=3meegeven aanAnova()
Meer uitleg over effectcodering en type III toetsen kan je vinden in deze syllabus vanaf sectie 5.3. De uitleg gebeurt in de context van klassieke lineaire regressie (dus niet multilevel modellen) maar de principes zijn identiek.
Voor de eerste stap moeten we het model opnieuw fitten, waarbij we
een extra argument contrasts meegeven.
model.ia.effect <- lmer(formula = BL_posttest ~ BL_pretest + Conditie*Dyslexie + (1|Klas),
data=lezen,
contrasts=list(Conditie=contr.sum, Dyslexie=contr.sum))
Vervolgens voeren we een modelvergelijking uit waarbij we het model
met het interactie-effect vergelijken met het model zonder. We kiezen
steeds voor type=3.
library(car)
Anova(model.ia.effect, type=3)Analysis of Deviance Table (Type III Wald chisquare tests)
Response: BL_posttest
Chisq Df Pr(>Chisq)
(Intercept) 0.0578 1 0.8099
BL_pretest 132.1697 1 <2e-16 ***
Conditie 0.1735 1 0.6770
Dyslexie 0.2971 1 0.5857
Conditie:Dyslexie 0.2268 1 0.6339
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
De p-waarde is ongeveer gelijk aan 0.634. Dit is groter dan 0.05 waardoor we de nulhypothese niet kunnen verwerpen op het 5% significantieniveau.
Visualisatie
Het interactie-effect tussen Conditie en
Dyslexie kan eenvoudig gevisualiseerd worden met het
package effects. Je geeft Conditie*Dyslexie en
het model model.ia aan de functie effect() en
geeft dit op zijn beurt aan de functie plot().
install.packages('effects') # eenmalig
library(effects)
plot(effect('Conditie*Dyslexie', model.ia)) 
Omdat Conditie*Dyslexie een interactie-effect voorstelt
krijgen we meerdere plots terug. We zien in de output een plot van het
effect van Conditie op BL_posttest wanneer
Dyslexie de waarde Ja aanneemt en een plot van
het effect wanneer Dyslexie de waarde Nee
aanneemt. In lijn met de output van summary(model.ia)
stellen we vast dat het effect van Conditie iets groter
lijkt te zijn bij de leerlingen die Dyslexie hebben (al
werd er geen statistisch significant verschil tussen de effecten
vastgesteld).
6 Varianten en uitbreidingen
Het voorbeeld op deze pagina reikt handvaten aan om zelf multilevel modellen te fitten. Er zijn echter veel concrete situaties denkbaar waarin je de hier gedemonstreerde modellen en bijhorende code zal moeten aanpassen zodat ze jouw concrete onderzoeksvraag kan helpen beantwoorden.
Opties om de hierboven gedemonstreerde technieken toepasbaar te maken op jouw specifieke casus zijn onder meer de volgende.
- Meer predictoren. Er is geen harde grens aan het aantal predictoren in een multilevel model. Het is dus mogelijk om (nog) meer predictoren toe te voegen aan het model. Het kan daarbij gaan om hoofdeffecten of interactie-effecten.
-
Meer niveaus. In het voorbeeld tot hiertoe waren de data
geclusterd op een hoger niveau (niveau 2, i.c. het klasniveau). Het is
ook mogelijk dat data nog verder geclusterd zijn op een nog hoger, derde
niveau.13 De klassen maken op hun beurt deel uit
van scholen. Als onderzoeker kan je vermoeden dat de testscores
(
BL_posttest) binnen scholen ook wat meer op elkaar lijken in vergelijking met alle leerlingen over alle scholen heen. -
Random slopes. In de demonstraties op deze pagina kon het intercept variëren per klas (of algemener geformuleerd, per cluster). In andere multilevel modellen is het ook mogelijk dat de grootte van een effect varieert per cluster. We spreken dan van een “random slope model”. Dergelijke modellen worden gefit wanneer niet voldaan is aan de assumpties van het random intercept model en/of wanneer er theoretische redenen zijn om te denken dat een effect verschillend is per cluster.
In random slope modellen hebben onderzoekers soms bovendien de bedoeling om de variabiliteit in de slopes tussen klassen te verklaren aan de hand van predictoren. Dan probeert men dus te verklaren waarom de slope in sommige clusters groter of kleiner is dan in andere clusters.
7 Referenties
Bogaert R., Merchie E., Van Ammel K. & Van Keer H. (2025). The impact of a tier 1 intervention on fifth and sixth graders’ lezen comprehension, lezen strategy use, and lezen motivation. Learning Disability Quarterly 48 (2), 102-115. https://doi.org/10.1177/07319487221145691
Chen D. & Chen J.K. (2021). Statistical regression modeling with R. Longitudinal and multi-level modeling. Springer. https://doi.org/10.1007/978-3-030-67583-7
Hox, J. (2010). Multilevel analysis: techniques and applications. 2nd ed. New York: Routledge.
Snijders T.A.B. & Bosker R.J. (2012). Multilevel analysis. An introduction to basic and advanced multilevel modeling. 2nd ed. London: Sage.
8 Voetnoten
Aangezien de interventie zelf op klasniveau gebeurt behoren alle leerlingen van dezelfde klas onvermijdelijk tot dezelfde conditie.↩︎
Leerlingen in de controlegroep legden ook twee testen af, net als leerlingen in de experimentele groep.↩︎
In klassieke lineaire regressie maken we de assumptie van onafhankelijke observaties. Wanneer de data geclusterd zijn is deze assumptie geschonden.↩︎
In sommige situaties is dat niet erg. Denk bijvoorbeeld aan PISA-studies, waar het niet altijd noodzakelijk is om te veralgemenen naar andere landen dan degene die opgenomen zijn.↩︎
De reden daarvoor is dat enkele functies in de beide packages dezelfde naam hebben. De functie van het laatst ingeladen package overschrijft de functie uit het voorgaande package, waardoor die laatste niet meer beschikbaar is. Hier gaat het om de functie
lmer().↩︎De ICC kan ook geïnterpreteerd worden als de correlatie van waarnemingen binnen de clusters. Waarom het begrip “correlatie” juist van toepassing is wordt uitgelegd in Chen & Chen (2021, p.31-32).↩︎
De visualisatie op deze pagina is gebouwd met het package
ggplot2.↩︎Het is perfect mogelijk om dit effect ook te laten variëren per klas (of algemener, op niveau 2). In dat geval spreken we van een random slope per klas.↩︎
In dit specifieke geval kan het resultaat van deze toets worden afgelezen uit de output van
summary()omdat er slechts 1 parameter is die codeert voor dit effect. Dit zal niet kunnen wanneer de categorische predictor meer dan 2 mogelijke waarden kan aannemen en er dus meer dan 1 parameter codeert voor het effect van die predictor. Bovendien is er in dit model geen interactie-effect aanwezig. Het toetsen van hoofdeffecten in modellen met een interactie-effect vereist extra aandacht en voorzichtigheid.↩︎Een equivalente manier om een interactie te coderen is met een
:tussenConditieenDyslexieen daarnaast explicietConditieenDyslexieals hoofdeffect toevoegen met+.↩︎Wanneer categorische variabelen in een model aanwezig zijn, is het altijd belangrijk om op de hoogte te zijn van de details van de codering, en in het bijzonder van het gekozen referentieniveau.↩︎
In dit specifieke geval kan het resultaat van deze toets worden afgelezen uit de output van
summary()omdat er slechts 1 parameter is die codeert voor dit effect. Dit zal niet kunnen wanneer de categorische predictoren in het interactie-effect meer dan 2 mogelijke waarden kunnen aannemen en er dus meer dan 1 parameter codeert voor het interactie-effect.↩︎Of zelfs op een vierde, vijfde,… niveau.↩︎